rm(list = ls())
set.seed(5168)

	library(foreign)
	library(mnormt)
	library(ordinal)
  library(MASS)

  data <- read.dta('Netherlands (2012).dta')

  data <- subset(data, data$Party != 1) 
  data <- subset(data, data$Party != 2) 
  data <- subset(data, data$Party != 4) 
  data <- subset(data, data$Party != 5) 
  data <- subset(data, data$Party != 6) 
  data <- subset(data, data$Party != 7) 
  data <- subset(data, data$Party != 8) 
  data <- subset(data, data$Party != 9) 
  data <- subset(data, data$Party != 10) 

  model <- polr(as.factor(Resp_attribution) ~ Prime_minister + Cabinet + Opposition + Perc_seats + Prime_minister * Perc_seats + Cabinet * Perc_seats + Opposition*Perc_seats, data = data, Hess=TRUE, method = 'probit')
	summary(model)

	pars <- as.matrix(model$coefficients)
  pars[8:11] <- as.matrix(model$zeta)     #since this is using a different package I need to add the intercepts to the vector that takes all of the coefficients
	vce <- vcov(model)


	simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

	summary(model)
	summary(simCoef)

	pmHi <- c()
	pmLo <- c()
	pmMed <- c()

  partnerHi <- c()
  partnerLo <- c()
  partnerMed <- c()

	supHi <- c()
	supLo <- c()
	supMed <- c()

	oppHi <- c()
	oppLo <- c()
	oppMed <- c()

  otherHi <- c()
  otherLo <- c()
  otherMed <- c()


	for(i in 0:50){
	## Prime minister 

		agg <- simCoef[, 1] * 1 + 
				simCoef[, 2] * 0 + 
				simCoef[, 3] * 0 + 
				simCoef[, 4] * i + 
				simCoef[, 5] * i + 
				simCoef[, 6] * 0 + 
				simCoef[, 7] * 0 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		pmHi[i + 1] <- quantile(pred, 0.975)
		pmLo[i + 1] <- quantile(pred, 0.025)
		pmMed[i + 1] <- quantile(pred, 0.5)


	## partner
	agg <- simCoef[, 1] * 0 + 
	  simCoef[, 2] * 1 + 
	  simCoef[, 3] * 0 + 
	  simCoef[, 4] * i + 
	  simCoef[, 5] * 0 + 
	  simCoef[, 6] * i + 
	  simCoef[, 7] * 0  
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		partnerHi[i + 1] <- quantile(pred, 0.975)
		partnerLo[i + 1] <- quantile(pred, 0.025)
		partnerMed[i + 1] <- quantile(pred, 0.5)

	## opposition
	agg <- simCoef[, 1] * 0 + 
	  simCoef[, 2] * 0 + 
	  simCoef[, 3] * 1 + 
	  simCoef[, 4] * i + 
	  simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * i 
		
		pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
		pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
		pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
		pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
		pred5 <- 1 - pred1 - pred2 - pred3 - pred4
		
		pred <- pred5 * 5 + 
			pred4 * 4 + 
			pred3 * 3 + 
			pred2 * 2 + 
			pred1

		oppHi[i + 1] <- quantile(pred, 0.975)
	  oppLo[i + 1] <- quantile(pred, 0.025)
	  oppMed[i + 1] <- quantile(pred, 0.5)

    ## No seats
	agg <- simCoef[, 1] * 0 + 
	  simCoef[, 2] * 0 + 
	  simCoef[, 3] * 0 + 
	  simCoef[, 4] * i + 
	  simCoef[, 5] * 0 + 
	  simCoef[, 6] * 0 + 
	  simCoef[, 7] * 0 
    
    pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
    pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
    pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
    pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
    pred5 <- 1 - pred1 - pred2 - pred3 - pred4
    
    pred <- pred5 * 5 + 
      pred4 * 4 + 
      pred3 * 3 + 
      pred2 * 2 + 
      pred1
    
    otherHi[i + 1] <- quantile(pred, 0.975)
    otherLo[i + 1] <- quantile(pred, 0.025)
    otherMed[i + 1] <- quantile(pred, 0.5)
    
    }


## now paint a pretty picture.... ##

	seats <- seq(0, 50, 1)        
  	
	#Now create the figure for Denmark
	par(mar=c(3.6,4,1,2))      
	
	plot(1, 1, type = 'n',
		xlim = c(0, 50),
		ylim = c(1, 5),
		xlab = '',
		ylab = '',
		main = '',
		axes = FALSE)
	axis(1,at=seq(0,50,25),labels=T)      
	axis(2)
	
	lines(seats, partnerMed, col = rgb(0, 0, 1, 0.2), lwd = 3)

	seatHi <- quantile(data$Perc_seats[data$Cabinet == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$Perc_seats[data$Cabinet == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], partnerMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 1, 0.3), lwd = 15)

	lines(seats, oppMed, col = rgb(0, 0, 0, 0.2), lwd = 3)

	seatHi <- quantile(data$Perc_seats[data$Opposition == 1], 0.9, na.rm = TRUE)
	seatLo <- quantile(data$Perc_seats[data$Opposition == 1], 0.1, na.rm = TRUE)
	lines(seats[seats >= seatLo & seats <= seatHi], oppMed[seats >= seatLo & seats <= seatHi], col = rgb(0, 0, 0, 0.3), lwd = 15)

	mtext("Responsibility attribution", side=2, line=2.5, cex=1.5)         
	mtext("Perceived seat %", side=1, line=2.2, cex=1.5)                 


############
###LEGEND###
############

legend('bottomright', legend = c('Partner', 'Opposition'), fill = c(rgb(0, 0, 1, 0.3), rgb(0, 0, 0, 0.3)), bty = 'n', border = NA, , cex=1.25)

###################
###EXPORT FIGURE###
###################

dev.copy(png,'Figure 7b NL 2012 - Effect of support party categorization.png', height = 490, width = 633)	
dev.off()  
